!> author: 左志华
!> date: 2022-06-24
!>
!> 气体状态方程
module sph_eos

    use sph_kinds, only: rk
    use sph_math, only: g

    procedure(p_water_fcn), pointer :: p_water => p_water_tait
    real(rk) :: c0                  !! 人工声速
    integer :: eos_type             !! 液体状态方程类型 <br>
                                    !! 1. tait
                                    !! 2. monaghan
                                    !! 3. morris
    real(rk) :: B                   !! 弱可压 EOS 参数, 默认值 ?
    real(rk) :: H_SWL               !!

    real(rk), parameter :: rho0 = 1000.0_rk
    integer, parameter :: gamma = 7

    abstract interface

        !> EoS for water <br>
        !> 水的状态方程
        pure subroutine p_water_fcn(rho, p, c)
            import :: rk
            real(rk), intent(in) :: rho     !! density <br>
                                            !! 密度
            real(rk), intent(out) :: p      !! pressure <br>
                                            !! 压力
            real(rk), intent(out) :: c      !! sound speed <br>
                                            !! 声速
        end subroutine p_water_fcn

    end interface

contains

    !> 应用状态方程通过密度和能量计算理想气体压力和声速
    pure subroutine p_ideal_gas(rho, u, p, c)
        real(rk), intent(in) :: rho     !! 密度
        real(rk), intent(in) :: u       !! 内部能量
        real(rk), intent(out) :: p      !! 压力
        real(rk), intent(out) :: c      !! 声速
        real(rk), parameter :: gamma = 1.4_rk

        associate (tmp => (gamma - 1.0_rk) * u)
            p = rho * tmp
            c = sqrt(gamma * tmp)
        end associate

    end subroutine p_ideal_gas

    !> Suitable for free surface flow <br>
    !> Tait simplified equation of state, 适用于自由面流动 <br>
    !> 孙鹏楠, 2018, 博士毕业论文, equ(2-12)
    pure subroutine p_water_tait(rho, p, c)
        real(rk), intent(in) :: rho     !! density <br>
                                        !! 密度
        real(rk), intent(out) :: p      !! pressure <br>
                                        !! 压力
        real(rk), intent(out) :: c      !! sound speed <br>
                                        !! 声速

        c = c0
        p = c*c*(rho - rho0)

    end subroutine p_water_tait

    !> Not suitable for free surface flow <br>
    !> 不适用于自由面流动
    !> @note 密度决定液体压强
    pure subroutine p_water_morris(rho, p, c)
        real(rk), intent(in) :: rho     !! density <br>
                                        !! 密度
        real(rk), intent(out) :: p      !! pressure <br>
                                        !! 压力
        real(rk), intent(out) :: c      !! sound speed <br>
                                        !! 声速

        c = c0
        p = c*c*rho

    end subroutine p_water_morris

    !> Suitable for free surface flow <br>
    !> 自由面流动
    pure subroutine p_water_monaghan(rho, p, c)
        real(rk), intent(in) :: rho     !! density <br>
                                        !! 密度
        real(rk), intent(out) :: p      !! pressure <br>
                                        !! 压力
        real(rk), intent(out) :: c      !! sound speed <br>
                                        !! 声速

        c = sqrt(B*gamma/rho0)*(rho/rho0)**3
        p = B*((rho/rho0)**gamma - 1)

    end subroutine p_water_monaghan

    !> Calculate sound speed <br>
    !> 计算声速
    pure subroutine get_sound_speed(h, c)
        real(rk), intent(in) :: h       !! specific enthalpy <br>
                                        !! 特定热量
        real(rk), intent(out) :: c      !! sound speed <br>
                                        !! 声速

        c = 10*sqrt(2*g*h)

    end subroutine get_sound_speed

end module sph_eos
